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ABSTRACT 

We analyze the linear stability of a stalled accretion shock in a perfect gas with 



O ! a parametrized cooling function L oc p^~"P°'. The instability is dominated by the 

^ \ 1 = 1 mode if the shock radius exceeds 2 — 3 times the accretor radius, depending 

"Th I on the parameters of the cooling function. The growth rate and oscillation period 

are comparable to those observed in the numerical simulations of Blondin & 
O ; Mezzacappa (2006). 

^ \ The instability mechanism is analyzed by separately measuring the efficiencies of 

' the purely acoustic cycle and the advective-acoustic cycle. These efficiencies are 

estimated directly from the eigenspectrum, and also through a WKB analysis 
^ ' in the high frequency limit. Both methods prove that the advective-acoustic 

03 I cycle is unstable, and that the purely acoustic cycle is stable. Extrapolating 

these results to low frequency leads us to interpret the dominant mode as an 
advective-acoustic instability, different from the purely acoustic interpretation of 
Blondin & Mezzacappa (2006). 

A simplified characterization of the instability is proposed, based on an advective- 
acoustic cycle between the shock and the radius ry where the velocity gradients 
of the stationary flow are strongest. The importance of the coupling region in 
this mechanism calls for a better understanding of the conditions for an efficient 
advective-acoustic coupling in a decelerated, nonadiabatic flow, in order to extend 
these results to core-collapse supernovae. 



Subject headings: accretion - hydrodynamics - instabilities - shock waves - su- 
pernovae 
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1. Introduction 

The recent discovery of a strong / = 1 instability of stalled accretion shocks in the 
context of core collapse supernovae (Blondin et al. 2003, Scheck et al. 2004, Ohnishi et al. 
2006, Burrows et al. 2006) has revived the interest in the fundamental stability properties 
of accretion shocks. This instability could be a major ingredient in the mechanism of ac- 
celeration of neutron stars (Scheck et al. 2004, Janka et al. 2004, Scheck et al. 2006a). 
It was also considered as a means to instigate g-mode dipole oscillations of the accreting 
neutron star (Burrows et al. 2006). While most of these authors recognized the presence of 
an advective-acoustic cycle similar to the one found by Foglizzo (2002, hereafter F02) in a 
different context, Blondin & Mezzacappa (2006, hereafter BM06) challenged this interpreta- 
tion and advocated a purely acoustic mechanism. Understanding the mechanism at work in 
this instability is a crucial step towards correctly extrapolating its consequences in a more 
realistic astrophysical situation. 

The physics of the advective-acoustic cycle is based on the linear coupling between 
acoustic and advectcd perturbations through the flow gradients: both entropy and vorticity 
perturbations act as source terms for the acoustic wave equation (Foglizzo 2001, hereafter 
FOl, and Foglizzo & Galletti 2003). This has been known for decades in the field of jet 
engines, since the pioneering works of Candel (1972), Howe (1975), Marble & Candel (1977) 
and Abouseif et al. (1984). In the subsonic flow below the stalled shock, this linear couphng 
is due to the gradients associated with the convergence of the flow, its deceleration, gravity 
and cooling. The interaction between the shock and the flow gradients gives birth to an 
advective-acoustic cycle, in which an advectcd perturbation generates a pressure feedback 
which triggers, at the shock, a new advected perturbation. Although Foglizzo (2002) proved 
the instability of the advective-acoustic cycle in an accelerated isothermal flow, the fate of 
such cycles in a cooled decelerated flow is still an open question: can they account for the 
instability observed in the simulations of BM06 ? 

The first aim of the present study is to clarify the instability mechanism at work, using 
perturbative techniques. The accretion flow is idealized as a perfect gas passing through 
a stationary shock, and subject to cooling processes schematically described by a cooling 
function L oc fJ^~°^P°', mimicking in the simplest manner the neutrino cooling in the core- 
collapse context. It is the flrst time that a linear approach has been used to understand the 
mechanism of this nonradial instability in the core-collapse context. 

A flrst step is to conflrm that the dominating I — 1 mode identifled by BM06 in the hnear 
phase of their numerical simulation indeed corresponds to the most unstable eigenmode of 
the linear problem. 
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Beyond the determination of the eigenspectrum and the vahdation of numerical simula- 
tions, wc wish to address the question of the instability mechanism, using techniques similar 
to F02. These techniques allow for a direct interpretation of the full eigenspectrum in terms 
of the efficiencies of the acoustic cycle and the advcctive-acoustic cycle. Alternatively, these 
efficiencies can also be computed in the WKB approximation. We wish to use both methods 
in order to check whether the instability is of acoustic or advective-acoustic nature. Under- 
standing the nature of the instability leads us to construct, in a companion paper (Foglizzo 
et al. 2006), a simple toy model which can be solved analytically and allows us to reach a 
fundamental understanding of some of the properties of the instability. 

The present paper is organized as follows: the boundary value problem associated with 
this stalled accretion shock is described in Sect. 2, where we estabhsh the boundary conditions 
at the shock and compare them to those used by Houck & Chevaher (1992, hereafter HC92) 
in the different context of supernova fallback. We determine in Sect. 3 the eigenfrequencies 
of the flows studied by BM06, and compare them with the linear phase of their numerical 
simulations. Then we investigate in Sect. 4 the mechanism responsible for this instability 
using the techniques of F02. The purely acoustic cycle is shown to be stable, and the 
advective-acoustic cycle is shown to be unstable with respect to / = 1 perturbations, in the 
range of validity of our approximations. These results are extrapolated to very low frequency 
perturbations in Sect. 5. The arguments of BM06 are reconciled with the advective-acoustic 
interpretation of the instability in Sect. 6. Results are summarized in Sect. 7. 



2. Formulation of the eigenvalue problem 

2.1. Description of the stationary flow 

We consider the radial accretion of a perfect gas with an adiabatic index 7 = 4/3, 
decelerated through a stationary shock at a radius Tgh, accreting on the hard surface of a 
neutron star of mass M and radius r*. The self-gravity of the accreting gas is neglected. 
The cooling function C is defined as a parametrized function of density p and pressure P as 
in HC92: 

C oc p^-"P", (1) 

which allows us to mimic the effect of neutrino cooling using the same prescriptions a = 3/2, 
/3 = 5/2orQ; = 6, /3 = las BM06. Neutrino heating, and the associated effect of convection 
is ignored in the present study. 

The equation of continuity, the Euler equation and the entropy equation defining the sta- 
tionary flow between the shock and the accretor are written in spherical coordinates as 
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Fig. 1. — Comparison of the stationary flows associated with the coohng parameters a = 6, 
(3 = 1 (full lines) and a = 3/2, f3 = 5/2 (dashed lines). The velocity and entropy profiles 
(upper plot) are similar in the outer parts of these two flows. The advection time Tadv from 
the shock to the accretor surface is flnite ii a < (3, and inflnite ii a > f3. The velocity 
gradient reaches a maximum at some intermediate height noted ry > if a > /3, whereas 
ry = if q; < p. The strength of cooling is chosen such that ry/rgh = 0.2 in both flows. 
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follows: 

|:(r» = 0, (2) 

d fv^ GM\ C 

or \2 7 — 1 r J pv 

- = -. (4) 

where S* is a dimensionless measure of the entropy defined as the following function of 
pressure and density, normalized by their values Pgh, psh immediately after the shock: 



7 



^ log 



P_ (Psh 
Ps^ V P 



(5) 



Note that the pressure force in the Euler equation (3) has been transformed using both this 
definition of S, and the sound speed c defined by = 7P/ p: 

VP / \ c2 

--vs. (6) 



p V7-1/ 7 

The shock is assumed to be adiabatic. Following HC92 and BM06, we assume that the 
pre-shock velocity Vi of the incoming gas is close to free fall: Vi ^ = —{2GM/rshY^'^, and 
that the gas is cold: A^i ^ 1. The Mach number J\4 = l^l/c is defined as a positive number. 

The assumption of stationarity required by the linear approach introduces a mathe- 
matical singularity at the surface of the accretor, where w(r*) = 0: the density diverges 
according to Eq. (2) and the sound speed decreases to zero. Such pathologies are common 
in linear studies of cooled accretion on a hard surface (e.g. from Chevalier & Imamura 1982 
to Saxton 2002), whatever the cooling function. For the cooling function considered here, 
two regimes can be distinguished depending on the sign of a — /3: 

(i) If a — j3 < 0, the cooling efficiency increases as the gas cools down, leading to a 
cooling runaway. The potential energy is negligible compared to the cooling losses, and the 
advection time to the accretor surface is finite. 

(ii) If a — /9 > 0, the cooling efficiency decreases as the gas cools down. The potential 
energy is comparable to the cooling losses, and the gas takes an infinite time to reach the 
surface. 

These two regimes are illustrated in Fig. 1, for the two set of cooling parameters used 
by BM06. The lower plot shows the velocity gradient of the stationary flow, which can 
participate to couple vorticity and acoustic perturbations (F02, Foghzzo & Galletti 2003). 
The velocity gradient reaches a maximum at a radius noted ry- Note that if a < /3, this 
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maximum is reached on the accretor surface (ry — r*). The strength of coohng in Fig. 1 was 
chosen such that the rgh/Tv = 5 in both flows. As noticed by BM06, the two flows are very 
similar in their outer parts. 

We find it convenient to use the variable logTW rather than the radius r in order to solve 
numerically the differential system in the cooling layer near the accretor surface. Integration 
is stopped just above the accretor surface, when the Mach number reaches 10~^. 



2.2. Differential system ruling the perturbed flow 

The flow is perturbed in 3-D using spherical coordinates. The complex frequency u = 
{uJr,uJi) of the perturbations is defined such that its real part ujr defines the oscillation 
frequency, and its imaginary part uji defines the growth rate. The perturbation of velocity 
6vr, Sve, 6v:f, density 6p, sound speed 6c, and entropy 6S are used to define new perturbative 
functions /, h, SK, which enable a compact formulation of the differential system once 
projected on spherical harmonics Yj^{9, (p): 



f = v5vr + 

5Vr 



7-1 
2 5c 



c5c, 



9 



V 7 — 1 c 



5K = r v.W X 5w + 



7 



(7) 
(8) 
(9) 



where 5w = V x Sv is the perturbation of vorticity. The resulting differential system is 
independent of the azimuthal number m: 
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(11) 
(12) 

(13) 
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where the quantity //^ used in Eq. (11) is defined by: 

1^' = 1- ^ ^ ; 1-M^ ■ 14 

The lengthy equations describing how the functions /, h, 6S, 6K translate into the classical 
quantities 5v, 6p, 6P, and the explicit expression of S{C/ pv) and 5{C/ Pv) in terms of /, /i, 5S, 
are written in Appendix A. 



2.3. Boundary conditions at the shock 



The boundary conditions arc established by writing the conservation laws in the frame 
of the perturbed shock. The derivation of these boundary conditions is shown in detail in 
Appendix B, especially since we do not use the special system of spatial coordinates used by 
HC92. These moving spatial coordinates were introduced in perturbative studies of accretion 
shocks on white dwarfs (Chevaher & Imamura 1982), such that the shock coordinate remains 
fixed even when the flow is perturbed. By contrast, our boundary conditions arc expressed 
at the radius rsh of the unperturbed shock, as functions of the displacement of the shock, 
associated with its velocity Aw = —iuj^C,. If the shock is strong {M-l-^ — {l ~ l)/27) and 
the incoming gas is in free fall {yi = vg): 
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(15) 
(16) 
(17) 
(18) 



These conditions can be translated into the classical quantities 6vr, Sp, 6P: 
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(21) 



The transverse components Sve, Sv^ of the velocity perturbation after the shock are related 
to the nonspherical deformation of the shock through: 



Tsh sin 9 d(fi 



(22) 
(23) 



The divergence of the transverse velocity can be projected on the spherical harmonics YJ^ 
and satisfies the following boundary condition for a strong shock: 



sin 9 



d_ 
89 



sva.95v0 



9 , ■ 



2/(/ + 1) 
■ 7-1 



(24) 



Equations (19), (20) and (21) arc in perfect agreement with Eqs. (50), (52) and (53) of HC92, 
which can be recovered by taking into account the gradients of the stationary flow quantities 
V, p and P, using a Taylor expansion between Tgh and Tgh + A^. 

Transverse velocities at the shock are precluded by Eq. (51) of HC92, in contradiction with 
our Eqs. (22-23) and (24), and also with other linear studies involving transverse perturba- 
tions such as Bertshinger (1986), Imamura et al. (1996), or Saxton & Wu (1999). 
As a consequence, the stability results reported by HC92 for nonradial perturbations should 
be considered as questionable, while their results concerning radial perturbations should be 
unaffected. 



3. Numerical determination of the eigenfrequencies 

3.1. Numerical method 

The differential system is solved by integrating over the variable logM. from the shock 
down to the accretor surface, at a point where M. = 10^^ if the advection time is finite 
{a < P). As an illustration, the radial shape of the eigenfunctions associated with the most 
unstable I — 1 mode in the flow with a — 3/2, (3 — 5/2, Tgh/r* = 10, is shown in the upper 
plot of Fig. 2. Using log as a variable allows us to compute the eigenfunctions down to 
the singular accretor surface. Structures are visible down to ~ 10~^, corresponding to 
radial scales which are much too small to be accessible to existing numerical simulations 



-9- 




1# 0,0001 0^001 M 0,01 m 



Fig. 2. — In the upper plot, radial profiles of the perturbations of \6K\, \SS\ and \6P\/P for 
the most unstable / = 1 mode, in the fiow with a = 3/2, f3 = 5/2, Vsh/r^ = 10. The radial 
profiles are shown as a function of the Mach number Ai, with the corresponding value of the 
fractional radius shown on the right axis. The lower plot shows the influence of the radius rbc, 
and associated Mach number At bo of the lower boundary, on the complex eigenfrequency 
Lo: the very fine structure of the perturbations for AI < 4 x lO""^ can be safely neglected. 



-10- 



(< 10~^rsh)- Luckily, these scales do not need to be resolved to measure the correct eigenfre- 
quency, as shown in the lower plot of Fig. 2. Varying the depth rtc at which the boundary 
condition Sv{ri,c) = is applied indicates that the regions of the flow where Ai < 4 x 10~^ 
(i.e. a fraction < 0.1% of the shock distance) have a negligible effect on the eigenfrequency. 
If the advcction time is infinite {a > f3), the boundary condition 6v/v = is applied at 
a radius where the advection time from the shock reaches 10 times the reference timescale 
('"sh — '"*)/|^sh|- We expect eigenmodes with a growth rate comparable to jt'shl/l^'sh — ^"*) to 
be insensitive to any ad vective- acoustic artifact associated with this numerical prescription, 
because it occurs on a longer timescale than the instability. This expectation is validated by 
checking that the eigenfrequencies are unchanged by increasing the advection depth. 
Once an eigenfrequency is found for a given intensity of the cooling function, it is tracked 
in the complex plane using the Newton-Raphson method. Eigenfrequencies are expressed in 
units of |t'sh|/(?^sh — '>^*) throughout the paper. 

3.2. Compcirison with the eigenfrequencies estimated by BM06 

BM06 validated their numerical code in 1-D by comparing the eigenfrequency measured 
in the linear stage to the eigenfrequencies determined by HC92 for the mode / = 0, with 
a = 5/2, [3 = 5/2. Our calculation confirms this validation, as shown in Fig. 3. Incidentally, 
Fig. 3 shows that this flow is much more unstable to / = 1 perturbations than to radial ones. 
The nonspherical axisymmetric calculations of BM06 have been performed in a flow in which 
a = 3/2, /? = 5/2, which was not considered by HC92. The eigenfrequencies of this flow are 
shown in Fig. 4, with a globally acceptable agreement. The agreement seems significantly 
better for the I — mode than nonradial ones. Fig. 4 also shows that higher harmonics 
dominate the instability if the shock is far enough. For Tgh/r* = 5, the instability of the 
mode / = 1 should be dominated by the first harmonics, whereas the mode I — 2 should be 
dominated by the second harmonics. 

Fig. 5 shows the expected growth rate in the flow with a = 6, i3 = 1 also considered by 
BM06, with a comparable agreement. Both Fig. 4 and 5 enables us to evaluate the accuracy of 
numerical simulations: both the oscillation period and the growth time should be considered 
with a typical 30% uncertainty. 

Some of the discrepancies may be attributed to difficulties in disentangling higher harmonics 
which have a similar growth rate. Whether part of this discrepancy could be attributed to 
numerical viscosity or not depends to some extent on our understanding of the instability 
mechanism. At low frequency, a purely acoustic mechanism should be barely sensitive to 
numerical viscosity. By contrast a mechanism involving the advection of vorticity waves 
towards regions of small velocity may be more sensitive to numerical viscosity. This could 
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Fig. 3. — Oscillation frequency ujr and growth rate Ui of the fundamental modes / = 0, 1, for 
a = (3 = 5/2. Also shown are the / = eigenfrequencies computed by HC92 (circles) and 
those measured in the 1-D simulation of BM06 (diamonds). The mode / = 1 is always more 
unstable than the radial mode in this flow. 
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Fig. 4. — Eigenfrequencies for a cooling law with a — 3/2 and (3 = 5/2, corresponding to 
the modes I — (dotted line), I — 1 (full hne) and I — 2 (dashed hne). The fundamental 
mode is plotted with thick black lines. Harmonics are shown with thiner grey lines. The 
eigenfrequencies determined from BM06 are shown as triangles {I — 0), squares {I — 1) and 
circles (1 — 2). 
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Fig. 5. — Same as Fig. 4, for the cooling law a — 6 and (3 — 1. 
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Fig. 6. — Eigenfrequencies corresponding to the fundamental modes < Z < 40, for a = 3/2 
and (3 = 5/2. This plot indicates that the most unstable mode is always nonradial. The 
larger the shock radius, the smaller the degree Z > 1 of the most unstable mode. 
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contribute to explain why the growth rate measured in the nonradial numerical simulations 
seems systematically lower than that predicted by the linear analysis. The radial profiles of 
SS and SK, displayed in Fig. 2, show structures in the region of ~ 10~^ which cover a 
fraction ~ 2% of the shock distance (see also the lower plot of Fig. 18). Given the difficulty 
of advecting vorticity waves in a grid based code, a grid size of 0.1 — 0.2% of the shock 
distance might be desirable near the accretor surface. Whether the 300-450 radial zones 
used in the simulations of BM06 are sufficient or not could be easily checked by performing 
new numerical simulations on a finer grid. 

BM06 have pointed out the resemblance between the instability in these two flows 
despite the different cooling function. Although we agree on the resemblance between these 
two flows when the shock radius is large, Fig. 5 reveals some signiflcant differences for smaller 
shock radii: 

- the mode / = 1 is always the most unstable mode if a = 6, /3 = 1. By contrast, if 
a — 3/2, /3 = 5/2, the instability may be dominated by perturbations with a higher degree 
Z > 2 if the shock radius is smaller than 2r*, 

- the flow with a = 6, = 1 is stable if the shock radius is shorter than ~ 2.5r*, whereas 
the flow with a = 3/2, /? = 5/2 is unstable whatever the shock distance. 

A more detailed investigation of the instabihty at small shock distance when a — 3/2, 
P — 5/2, illustrated by Fig. 6, suggests that the most unstable mode corresponds to an 
azimuthal structure with a size comparable to the shock distance. 

For both coohng functions, the dimensionless growth rate of the instabihty is at best of the 
order of 0.7 - 0.8: 

< — - — 

'"sh ^* 

This similar growth rate can be viewed as a hint of a common physical mechanism of insta- 
bility, which we investigate in the next section. 

4. Determination of the instability mechanism 

4.1. Presence of oscillations in the eigenspectrum 

The oscillation period and growth time associated with the most unstable eigenmode, 
as determined in Figs. 4 and 5, should be a signature of the instability mechanism. Unfor- 
tunately, our understanding of the possible instabilities is not deep enough to allow for a 
direct and conclusive interpretation of these timescales. The linear stability analysis can be 
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Fig. 7. — Eigenfrequencies computed from the boundary value problem in a flow with a = 
3/2, P = 5/2 (upper plot) and a — 6, (3 — 1 (lower plot), for different shock radii rgh/rv, for 
I — 1 perturbations. 
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helpful in determining the underlying mechanism, using the many other eigenfrequencies of 
the eigenspectrum. A global view of the eigenspectrum associated with / = 1 perturbations 
is shown in Fig. 7 for both cooling functions, with different shock distances. Figure 7 sug- 
gests that the larger the shock radius, the more numerous the unstable modes. A striking 
feature of these eigenspectra is the oscillations of the growth rate, which are easier to identify 
when the number of well defined eigenmodes is large, i.e. when the shock radius is large. 
According to Fig. 7 , this identification requires rgh/^^v > iO ii a — 6, (3 — 1, whereas the 
oscillations are already visible for Tgh/r* = 5 if a = 3/2, P — 5/2. These eigenspectrum 
oscillations are reminiscent of the oscillations visible in the eigenspectrum of the isothermal 
accretion fiow accelerated towards a black hole (Fig. 4 of F02). They were explained by F02 
as the consequence of the influence of a purely acoustic cycle interacting either constructively 
or destructively with the advective-acoustic cycle. The efficiencies of these two cycles are 
measured in the next Section, after recalling the formalism associated with these cycles. 

4.2. Calculation of the efficiencies Q, TZ of the advective-acoustic and purely 

acoustic cycles 

4-2.1. Advective-acoustic and purely acoustic cycles 

F02 developed a formalism in order to describe the advective-acoustic cycle, stable or 
not, below a stationary shock in a radial accretion fiow onto a black- hole. The same for- 
malism can be applied to a decelerated accretion fiow onto a hard surface. Let us recall 
that vorticity and entropy perturbations are advected at the velocity of the fiow, whereas 
pressure perturbations can propagate at the speed of sound. These two categories of per- 
turbations would be linearly independent in a uniform flow, but are coupled linearly if the 
flow is inhomogeneous. The interaction between the stationary shock and the fiow gradients 
gives birth to two cycles: 

- an advective-acoustic cycle, whose duration is noted tq. An advected perturbation 
of frequency ujr generates a pressure feedback which triggers, at the shock, a new advected 
perturbation, whose amplitude has changed by a factor Q{u!r) after one advective-acoustic 
cycle. 

- a purely acoustic cycle , whose duration is noted t-jz- An acoustic perturbation of 
frequency Ur propagating downward (not necessarily radially) produces a reflected (or re- 
fracted) perturbation reaching the shock and triggers a new pressure perturbation, whose 
amplitude has changed by a factor TZ{ujj) after one acoustic cycle. 

The "vortical-acoustic" instability studied by F02 is fundamentally nonradial because vor- 
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ticity perturbations are the only advected perturbations in an isothermal flow. In a gas 
with 7 = 4/3, the advective-acoustic cycles exist for both radial and non radial pertur- 
bations. A radial advective-acoustic cycle relies entirely on entropy perturbations. In an 
adiabatic flow, the acoustic feedback depends on the global increase of enthalpy in the post- 
shock flow (Foglizzo & Tagger 2000, FOl) and can be efficient enough to destabilize the 
Bondi-Hoyle-Lyttleton accretion (Foglizzo, Galletti & Ruffert 2005). The instability of this 
"entropic-acoustic" cycle (i.e. the mode Z = 0) is disfavoured in the core-collapse context 
because neutrino cooling precludes a strong adiabatic heating. 



4-2.2. Eigenfrequencies associated with the cycles 

The simplest formulation of the advective-acoustic instability corresponds to a situation 
where the purely acoustic cycle is negligible. The instability threshold then corresponds to 
|Q| = 1 and the growth rate cui can be approximated by 

u;, ~— log|Q|. (26) 

More generally, Foglizzo & Tagger (2000) showed that the purely acoustic cycle is not nec- 
essarily negligible and modifies Eq. (26) as follows: 

Q^i-^TQ ^ ^gia;r^ ^ ^ ^27) 

This equation describing the simultaneous existence of two cycles (QjTq) and {TZ,ttz) is 
symmetric: it can account for an advective-acoustic instability (|Q| > 1) as well as an 
hypothetical acoustic instability (if |7^| > 1). In the isothermal fiow studied by F02, the 
acoustic cycle is "weak" in the sense that the parameter e < 1: 



I7^l 



IS. 



I She 



< 1, (28) 



Assuming e < 1, F02 showed from Eq. (27) that the effect of the acoustic cycle is to either 
increase or decrease the growth raye cUi in the following range: 

l,og-^<..<l,ogJ^. (29) 

Tq 1 + 6 Tq 1 - e 

The case e > 1 would be exactly symmetric, by exchanging (Q, tq) and {TZ, t-ji) in Eqs. (28) 
and (29). 
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Fig. 8. — Comparison between the coupling efficiencies |Q|(a;), |7^|(tf) of the mode / = 1, 
computed in the WKB approximation (full and dashed lines) and the ones deduced from 
the eigenspectrum (circles and diamonds) using Eqs. (32-33), in a flow with a — 3/2 and 
/3 = 5/2, for different shock radii Tgh/r*. The agreement validates the formahsm associated 
with the two cycles. The acoustic cycle is stable {\Tl\ < 1), and the ad vective- acoustic cycle 
can be unstable (|Q| ~ 4 — 6). 
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Fig. 9. — Same as Fig. 8 with a different cooling function: a = 6 and (3 — 1. The efficiencies 
\Q\ and \7l\ are similar to those obtained with a = 3/2 and /3 = 5/2. 
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4-2.3. How to extract the cycles information directly from the eigenspectrum 

Assuming that the underlying mechanism is due to a superposition of cycles described 
by Eq. (27), an estimate of Q, TZ, tq, tji can be extracted directly from the oscillations 
observed in the eigenspectrum, in Fig. 7. This method enables us to identify two cycles, one 
slow and one fast, and measure their efficiencies. For the sake of simplicity, we choose to 
denote the fast cycle with the letter TZ, and the slow one with the letter Q (i.e. tq > Tn)- 
The identification of the slow cycle with the ad vective- acoustic mechanism, and the fast cycle 
with the purely acoustic mechanism will become unambiguous in Sect. 4.2.4, which will also 
validate the assumption of an Eq. (27) underlying the eigenspectrum. 

According to F02, the timescale tq of the slowest cycle is related to the frequency difference 
ujr{i + 1) — ujr{i) between two consecutive eigenmodes, whereas the shortest timescale t-ji 
is related to the frequency range Aa;^ of each oscillation. Denoting by rzosc the number of 
eigenmodes per oscillation, 

— ~ riosc- (31) 

For example, at first glance on the upper plot of Fig. 7, one can anticipate that Tq/t-ji ~ 5 
if Tsh/T* = 5, whereas tq/t-jz ~ 10 if rgh/T* = 30. 

According to Eqs. (54) and (55) of F02, the efficiencies |Q| and \TZ\ associated with the slow 
and fast cycles respectively, can be extracted directly from the eigenspectrum, by measuring 
the following parameters of the eigenspectrum oscillations: (i) the frequency range Aujr of 
each oscillation, (ii) the number rzosc of eigenmodes per oscillation, (iii) the amplitude Acui 
of the oscillations of the growth rate, and (iv) their average value cuf 

coshvr^ O- 

\Q\ = — TT it Ac. exp2nosc7r- — , (32) 

cosh(nosc - l)7r|^ AuJr 

sinhnoscTTx^ a;,- 
^ 2n-^. (33) 

cosh(nosc - l)7r|^ Ao;, 

The result of this method, applied to the / = 1 eigenmodes of Fig. 7, is shown in Fig. 8 for 
a < f3 and in Fig. 9 for a > (3. It shows that the fast cycle (diamonds) is always stable, 
whereas the slow cycle (circles) can be unstable. We now proceed to check the validity of 
these results by computing \ Q\ and \TZ\ using another method, which will establish that the 
fast stable cycle is purely acoustic, and the slow cycle is advective-acoustic. 
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4-2. 4- An alternate way to measure \ Q\ and \TZ\, in the WKB approximation 

Following F02, this method consists in determining the coupling efficiencies |Q| and 
1 7^ I as continuous functions of the perturbation frequency ujr- The efficiency of the acoustic 
cycle is decomposed into TZ = TZghR-s/, while the efficiency of the advective- acoustic cycle is 
decomposed into Q = QshQv, with TZsh, "T^v, Qsh and Qy being defined as follows: 

(i) when an acoustic wave of frequency Ur propagating outward reaches the shock, 
TZsh{i-^r) measures the efficiency of acoustic reflection, while Qshii-^r) measures the amount of 
advected perturbations (entropy/vorticity) produced by the shock. 

(ii) when an acoustic wave of frequency ujr propagates in the flow towards the accretor, 
7^v('^r) measures the amount of acoustic reflection. When an advected perturbation of fre- 
quency Ur is advected towards the accretor, Qv(co'r) measures the amount of acoustic waves 
propagating against the flow. Note that Qy and T^v are measured at a radius immediately 
below the shock, but do not involve the physics of the shock. 

This approach is an extension of the approach used by FOl, F02 in adiabatic and isothermal 
flows, where acoustic and advected perturbations are easily identified using a WKB approx- 
imation. The technicalities of the method is described in Appendix C and Appendix D. 
The WKB approach assumes that the wavelength of the advected and acoustic perturbation 
is shorter than the scale of the flow gradients just below the shock. This method is thus 
expected to be reliable at high frequency, and to break down at low frequency. 
Besides, the presence of cooling processes makes this method even more approximative, since 
we choose to neglect it in the immediate vicinity of the shock, for the sake of simplicity. Ne- 
glecting cooling is certainly justified when the shock is far enough from the accretor, as 
illustrated by the flat entropy proflle in the upper plot of Fig. 1, for rgh/ry = 5. 

The results of the WKB analysis appfied to the I — 1 mode are displayed as continuous 
and dashed lines in Figs. 8 and 9, together with the results of the flrst method. The agree- 
ment between the two methods is excellent at frequencies where |Q| and \TZ\ vary smoothly 
with frequency. Since the method based on the eigenspectrum requires several neighboring 
eigenmodes to determine \ Q\ and \TZ\, it is unable to correctly capture variations that are 
faster than the oscillation period Aa;^. This is particularly visible in the resulting estimate 
of \n\, in Fig. 8. 

The direct method used in Sect. 4.2.3 can be used to estimate the range of validity of the 
WKB approximation. A close inspection of Figs. 8 and 9 suggests that the WKB approxi- 
mation is acceptable for cUr > 5 in both flows. 
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Both methods indicate that the purely acoustic is stable, and is not even close to the 
instability threshold, with typical values \TZ\ < 0.5. If a purely acoustic mechanism were 
responsible for the low frequency instability, as proposed by BM06, this mechanism would 
have to be stable at higher frequency to be compatible with the results of our analysis. 

More importantly, our calculations prove that the advective-acoustic cycle is unstable 
in the range of frequencies accessible to our analysis, with an efficiency reaching \Q\ ~ 4 — 6. 
The next sections aim at characterizing the properties of this instability, by estimating 
the timescale tq of the cycle and the effective coupling radius Vcs (Sect. 4.3), and finding 
approximations for its growth rate (Sect. 4.4) and oscillation period (Sect. 4.5). Whether this 
instability mechanism may be responsible for the dominant low frequency mode is discussed 
in Sect. 5. 



4.3. Estimate of the timescale tq and effective coupling radius Ves of the 



The accurate determination of the cycle timescales tq and tti, using the velocity and 
sound speed profiles of the stationary flow, is not straightforward: the duration Tn of the 
acoustic cycle depends on the depth of the turning point, and is thus a function of both the 
frequency cUr and the order I. Similarly, cUr and I influence the depth at which an advected 
perturbation couples most efficiently to acoustic perturbations. The study of the advective- 
acoustic coupling in an adiabatic or isothermal flow showed that this coupling occurs all 
the way from the shock to the accretor (Eq. (25) of FOl, of Eq. (29) of F02). Its net effect 
viewed from the shock radius may be summarized as a feedback from a single effective radius 
of coupling Teff, associated to an advective-acoustic timescale tq. Let us define two radial 
reference timescales in the stationary flow: 

(i) the radial acoustic timescale Tg^{r) from the shock to the radius r and return: 



(ii) the radial advective-acoustic timescale raa(r) is defined as the advection timescale 
from the shock to the radius r, and acoustic return in the radial approximation: 



advective-acoustic cycle 
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Fig. 10. — Estimate of the cycle timescale tq, directly extracted from the eigenspcctrum 
(Eq. 30). It is measured in units of the radial advective-acoustic time ry = Taa('"v)) down to 
the radius ry of maximum velocity gradient. At low frequency, the cycle timescale is 20% 
shorter than Taaf^v)- 
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The non radial character of the acoustic feedback has a minor influence on raa(r), given the 
subsonic character of the flow. By contrast, BM06 noted that the acoustic time for a very 
nonradial, purely acoustic wave can be significantly longer than along a radial path: the 
shock circumference is indeed a factor tt longer than its diameter. 

In Fig. 10, we find it convenient to measure the cycle timescale tq deduced from Eq. (30), 
in units of ry = Taa('"v)! where ry is the radius defined in Sect. 2.1. where the velocity gra- 
dient is maximum. The globally good matching between these two timescales indicates that 
velocity gradients are an important ingredient for the advective-acoustic coupling responsible 
for the acoustic feedback, as in the vortical-acoustic instability studied in a isothermal con- 
text by F02. Note that temperature gradients may also contribute to the advective-acoustic 
coupling, as seen in the adiabatic study of FOl. 

We should keep in mind, however, that ry is defined as a local maximum, whereas our 
definition of Vcs is global in the sense that it involves both the local coupling efficiency and 
the advection/propagation of the perturbations between the shock and r^s- In consequence, 
the approximation ~ ry should be viewed as a guide to our intuition rather than a 
fundamental property. 

An advected perturbation of oscillation frequency Ur is most sensitive to fiow gradients 
whose lengthscale is shorter than the wavelength 27rv/uJr. Those associated with velocity 
scale like {dv/vdr)~^: 

ujr < 27r^. (36) 
dr 

According to Fig. 11, the velocity gradient in the fiow with a > j3 is very smooth and 
spread when the shock distance is short, whereas it gets sharper when the shock distance 
is large. This may explain, at least qualitatively, why the flow with a — (3 — 2 \s stable 
for Tsh/r* < 2.6 whereas the flow with Q; = 3/2,/3 = 5/2is unstable even for a small shock 
distance. 

Similarly, the divergence of the velocity gradient in the flow a < /3 is hkely to couple advected 
and acoustic perturbations at much higher frequencies than in the flow a > /?, where the 
velocity gradient is smoother. This qualitative argument may explain why the range of 
unstable frequencies is so much larger for a — 3/2, (3 = 5/2 than for a = 6, /3 = 1 (as visible 
in Fig. 7 and summarized in Fig. 16). 

Fig. 11 also illustrates the fact that velocity gradients are present all the way from the shock 
to ry, and may affect perturbations with a low frequency 0;^ < 5 according to Eq. (36). 
This threshold is comparable to our estimate of the threshold of the WKB approximation in 
Sect. 4.2.4. 

We interpret the 20% offset visible in Fig. 10 at low frequency as a consequence of the couphng 
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through the various flow gradients, including velocity, above ry- Keeping in mind this low 
frequency distortion, the timescale tv can be considered an acceptable approximation of tq 
for unstable / = 1 modes. 



4.4. First order approximation of the growth rate 

Finding a discrete set of eigenmodes is a laborious task, in comparison with the straight- 
forward calculation of |Q| and \TZ\ in the WKB approximation. In this respect, Eq. (29) 
provides us with a useful estimate of the growth rate if the parameters of the cycle tq, r-jz, 
\Q\ and \TZ\ are known. Although our estimate ry of tq is not fully satisfactory, we may 
continue along this direction in order to check what kind of accuracy can be reached. The 
effect of the acoustic cycle described by Eq. (29) depends through e on the ratio tq/t-ji 
(Eq. 28), which is measured by nose (Eq- 31). Despite the potentially significant difference 
between the timescales of radial and non radial acoustic modes, we choose to approximate r-ji 
by Tac(r*), and compare in Fig. 12 the number ^osc measured on the eigenspectra of Fig. 7 to 
the reference ratio Tv/Tac(r*). The global agreement is acceptable given the discreteness of 
the number of nodes, the gross approximation of ttz by Tac(r*), and the difficulty of identify- 
ing the oscillation, especially when Q varies rapidly near oOr ~ 90 and oUj. ~ 150 for a = 3/2, 
[3 = 5/2. The systematic decrease of riosc observed at low frequency is compatible with the 
decrease of tq/tv discussed in Sect. 4.3. 

The growth rate coi measured in the eigenspectra of Fig. 7 is then compared in Figs. 13 
and 14 to the value (log Q)/t^ (full line), and the expected range of influence of the acoustic 
cycle deduced from Eq. (29) (dotted lines). This comparison is interpreted as follows: 

(i) the amplitude of the eigenspectrum oscillation is very well matched by Eq. (29), 

(ii) the global shape of the eigenspectrum is globally very well reproduced 

(iii) as expected, some systematic discrepancies are observed at low frequency, concern- 
ing the flrst ~ 10 eigenmodes. This discrepancy can reach a factor 2 for the fundamental 
mode. 

We flrst conclude that the radius ry is an excellent approximation of the feedback radius at 
high frequency, for both cooling functions. This provides us with a rather simple description 
of the instability mechanism at work at high frequency. 

The discrepancy observed at very low frequency exceeds the 20% effect due to overestimat- 
ing Tq because the WKB approximation used to compute |Q| and \TZ\ ceases to be valid, as 
already pointed out in Sect. 4.2.4. Figs. 13 and 14 suggest that the WKB approximation 
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may only provide a gross estimate of the growth rate of the lowest frequency modes, within 
a factor 2. 



4.5. Oscillation timescale and efficiency \Q\ associated with the most unstable 

mode 

The relationship between the oscillation period of the fundamental mode and the timescale 
of the cycle is not obvious a priori, even if the acoustic cycle is neglected, because it depends 
on the phase 0q of the complex efficiency Q. Denoting by ujr{k) the frequency of the /c-th 
harmonic, the phase relation associated with Eq. (27) when the acoustic cycle is neglected 
leads to: 

UJr{k)TQ + (t)Q^2{k + l)T^. (37) 

A comparison of the frequency U0r{k = 0) of the fundamental mode with the frequency of 
the first harmonic ujr{k = 1), shown in Fig. 15, indicates that ujr{l) ~ 2uJr{0). We conclude 
that the phase of Q is negligible, and that the oscillation period of the fundamental mode is 
a good measure of the timescale tq at low frequency. 

The fundamental mode is not always the most unstable one among the / = 1 pertur- 
bations: Fig. 16 summarizes the frequency range of unstable modes for a = 3/2, /3 = 5/2 
(upper plot) and a — 6, (3 — 1 (lower plot). The most unstable mode may correspond to 
either the fundamental mode, the first or second harmonic as the shock radius increases 
from Tsh/rv = 1-5 to 100. For both sets of cooling parameters, the corresponding oscillation 
period 27r/u!r of the most unstable I — 1 mode would therefore be tq for 1.5 < rgh/ry < 5 — 6, 
tq/2 for 5 — 7 < Tsh/Tv < 33 — 39, and even tq/3 for Vgh/r^ > 33 — 39. 
An upper bound of the efficiency \Q\ associated with the low frequency modes may be 
estimated from Eq. (26) by neglecting the purely acoustic cycle and approximating tq ~ 
27r(A; + l)/ujr{k). The value of exp(cjjrQ) is shown in Fig. 17 for the first three eigenfre- 
quencies "f, "hi", and "h2", as a function of the shock radius. The actual value of Q at 
large shock radius is hkely to be intermediate between the curve "f ' and the curve "hi" in 
Fig. 17 for 10 < Tgh/rv < 30, as suggested by Fig. 7: indeed, the eigenspectrum oscillation 
at low frequency in Fig. 7 suggests that the first harmonics "hi" and "h2" benefit from a 
constructive infiuence of the acoustic cycle, whereas this influence seems destructive on the 
fundamental mode "f . In this respect, Fig. 17 indicates that |Q|wkb can be used as an 
acceptable guess of the efficiency \Q\ at low frequency. 

Besides, Fig. 17 suggests a slow increase of the efficiency \Q\ with the shock radius, with 
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remarkable similarity for the two sets of cooling parameters (full and dashed lines). 

As the shock radius is increased, the slow increase of the oscillation frequency of the most 

unstable mode (Fig. 16), and the increase of |Q| (Fig. 17) are not explained yet. 

5. Continuity argument for the advective-acoustic instability at low frequency 

Strictly speaking, we have demonstrated that the advective-acoustic cycle is responsible 
for the high frequency I — 1 instability of a stalled accretion when it is far enough from the 
accrctor, for two different types of cooling functions. What is the bearing of this demon- 
stration on the instability mechanism of the low frequency modes, when the shock radius 
is moderate ? Since the analysis of the fundamental / = 1 mode and its first harmonics is 
intractable through usual WKB techniques, one might argue that the instability mechanism 
of the fundamental mode is out of reach of the present study. The eigenspectrum of the 
flow with a — 6, (3 — 1 and Tgh/r^r = 5 (lower plot of Fig. 7) contains too few unstable 
modes to allow for the identification of oscillations, and the frequency of these eigenmodes is 
too low to allow for a WKB analysis. In this case, neither of the two methods described in 
Sect. 4.2 is able to compute \Q\ and \TZ\. Fig. 7, however, shows the continuity of the shape 
of the eigenspectrum, both with respect to frequency and with respect to the shock radius. 
The growth rate of the low frequency eigenmodes is only marginally larger than those at 
higher frequency. We feel there is no need to invoke a different instability mechanism at low 
frequency given the smooth distribution of growth rates. According to Fig. 7, decreasing 
the shock radius seems to reduce the cut-off frequency above which the modes are stable, 
but barely affects the growth rate of the low frequency modes. The continuity of the flow 
properties with respect to the shock radius is also apparent in the sequence of plots in the 
Figs. 8-9 and Figs. 13-14. 

Besides, the radial structure of the pressure perturbation looks very similar for all the 
harmonics displayed in Fig. 18, for two different shock radii. The eigenfunction of the 
low frequency, most unstable mode resembles those of higher frequency modes for which a 
reliable determination of the cycle efficiencies Q and TZ is possible. The similarity of the 
eigenfunctions suggests again a common instability mechanism for all these modes, namely 
the instability of the advective-acoustic cycle. 

Could the same kind of mechanism account for the slow instability of the fundamental 
I — mode, observed in Fig. 4 or 6 for a large shock radius Tgh/r* > 15 ? The eigenspectrum 
associated to the radial mode shows oscillations very similar to Fig. 7, except that all the 
higher harmonics are stable, even for a very large shock radius rgh/T* = 1000. These oscil- 
lations enable the direct computation of the efficiencies \Q\ and \TZ\ of the acoustic cycle. 
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The "entropic-acoustic" cycle is always stable |Q| < 1 at the high frequencies required by 
the WKB analysis, and the purely acoustic cycle is also stable \1Z\ < 1. Since the first 
harmonics seem always stable, the instability mechanism of the fundamental / = mode is 
more difficult to interpret than for I = 1 modes. We cannot exclude, however, that it may 
be due to a low frequency entropic-acoustic cycle destabilized by the temperature increase 
in the adiabatic part of the flow, when the shock radius is large. 



6. Discussion of the acoustic interpretation of BM06 

BM06 summarized their preference for a purely acoustic mechanism by the following 
three observations: 

(i) the advection time down to the accretor surface Tadv('"*) is larger than the oscillation 
period Tosc = 2Tr/ur, whereas a non radial path exists such that the acoustic timescale along 
this path coincides with Tosc (Fig. 6 of BM06), 

(ii) a standing I — 1 acoustic wave is clearly visible on the pressure profile. The particular 
radius r such that which could have been identified as a radius of effective 
ad vective- acoustic coupling, does not show any particular signature on the pressure profile, 

(iii) the instability seems independent of the flow near the accretor, given the resem- 
blance between the results with the two different cooling functions which differ only in the 

most inner regions. 

We believe that part of the confusion comes from the fact that for a post-shock Mach num- 
ber A4sh ~ 1/8^^"^ typical of a strong adiabatic shock with 7 = 4/3, the advection timescale 
^sh/^sh happens to be both longer than the radial acoustic timescale 2rsh/csh, and shorter 
than the surface acoustic time 27rrsh/csh: 

2rsh ^'sh 27rrsh /„ s 

^ < — < -. (38) 

Csh ^sh 

An advective-acoustic timescale of the order of rgh/^sh can thus be matched by an acoustic 
time along an ad-hoc nonradial path. This confusion between the advective and acoustic 
timescales could disappear if the postshock flow were very subsonic, with Aigh ^ (27r)~^ = 
0.16. Such conditions could be obtained in a flow involving a strong leak of energy at 
the shock, mimicking photodissociation in the same manner as Foglizzo, Scheck & Janka 
(2006). Unfortunately, preliminary calculations indicate that the shock distance is so much 
diminished by this energy loss that the I — 1 mode is stabilized. 



We argue below that the observations of BM06 do not contradict our description of the 
advective-acoustic instability. 
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- The resemblance between the two flows with different coohng functions (point iii) has 

been investigated in depth throughout this paper, showing that the instability is dominated 
in both flows by the advective-acoustic cycle occurring between the shock and the effective 
coupling radius r^s ~ ''"v- Some significant differences were also discussed, such as the smaller 
frequency cut-off in the fiow with a = 6, /? = 1 which could be due to the smoother profile 
of the velocity gradient (Fig. 11). 

- The existence of a standing pressure wave (point ii) cannot be held against the 
advective-acoustic cycle, as seen on Fig. 18. Most of the radial structure of the perturbed 
pressure shown in Fig. 18 is very localized in the first 10% of the shock distance, close to the 
accretor, and could be compatible with Fig. 7 of BM06 once projected on the / = 1 spherical 
harmonics. Let us recall that Teff is an "effective" coupling radius, as defined in Sect. 4.3, 
and certainly not the unique radius of acoustic emission. 

- As discussed in Sect. 4.3, our estimate of r^s ~ is not accurate enough to be a 
decisive test for our interpretation (point i). In a sense, our determination of r^s in Fig. 18 
could be considered as arbitrary as the determination of an "acoustic" path in Fig. 6 of 
BM06. Nevertheless, independently of estimating res, we did prove the instability of the 
advective-acoustic cycle at high frequency in Sect. 4, and showed the similarity with low 
frequency modes in Sect. 5. 

In view of the above detailed description of the advective-acoustic instabihty, together 
with the proven stability of the purely acoustic cycle at high frequency, the possibility that 
a purely acoustic mechanism may be responsible for the most unstable low frequency modes 
seems unlikely. 

7. Conclusions 

This work is the first characterization, through a linear study, of the advective-acoustic 
instability in a decelerated accretion flow involving cooling processes. Our formulation of the 
boundary conditions at the shock corrects an error in HC92 concerning non radial perturba- 
tions. The numerical solution of this problem confirms the existence of an / = 1 instability, 
as found in the numerical simulations of Blondin et al. (2003), for the two types of cooling 
functions studied by BM06. A detailed comparison of the growth time and the oscillation pe- 
riod of the dominant mode of the instability revealed discrepancies which can reach ~ 30%, 
possibly due to the numerical difficulty of advecting vorticity waves towards the region of 
deceleration without artificially damping them by numerical viscosity. The optimal grid size 
near the accretor surface is estimated as a fraction 0.1 — 0.2% of the shock distance from the 
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accretor. 

The main purpose of this study was to clarify the instabihty mechanism at work for I — 1 per- 
turbations, with the following results for both types of considered cooling functions a — 3/2, 
(5 = 5/2 and a = 6, = 1: 

(1) We have proven, for the first time, that an advective-acoustic instabihty of the I — 1 
mode takes place in a decelerated accretion flow involving cooling processes. The WKB 
approximation used in this proof required that the shock radius exceeds ~ 10 times the 
accretor radius. 

(2) The low frequency I — 1 instability occurring when the shock distance is moderate 
(rsh/Tv ^ 2) has also been interpreted as an advective-acoustic instability, using our con- 
clusion (1) together with a continuity argument: the instability of the low frequency modes 
can be interpreted in continuity with the instability at higher frequency in a series of flows 
with larger shock radii. This continuity argument is based on a comparison of both the 
eigenspectra (Fig. 7) and the eigenfunctions (Fig. 18) of the unstable modes. 

(3) The purely acoustic cycle is very stable {\R\ < 0.5) in the range of shock radii and 
frequencies allowed by our approximations. This result disfavours of the acoustic interpre- 
tation of BM06. 

(4) The efficiency of the advective-acoustic cycle is an increasing function of the shock 
distance, which reaches an efficiency |Q| ~ 4 — 5 for a shock radius rgh/Tv ^ 10. 

(5) We have proposed a simple approximation of the advective-acoustic cycle timescale 
tq based on ry, defined as the advection time from the shock to the radius ry where the 
velocity gradient is strongest, plus a radial acoustic feedback for the sake of simplicity. This 
estimate is quite accurate at high frequency, but overestimates the cycle timescale at low 
frequency by ~ 20%. 

(6) We have shown that the oscillation period of the fundamental mode is a measure 
of the timescale tq. The oscillation period of the most unstable mode is comparable to Tq, 
tq/2 or rQ/3 depending on the shock radius (Fig. 16). 

Our efforts to understand the instability mechanism at work in a simplified flow aim at 
guiding our intuition when interpreting more complex numerical simulations of astrophysical 
flows (i.e. Scheck et al. 2006b). The general description of the instability is globally satis- 
factory, but some fundamental questions are still unanswered. The following three questions 
may be answered by further studies: 

(i) What is the maximum efficiency |Q|max of an advective-acoustic cycle with cooling 
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? We observed in Figs. 17 that |Q|max increases slowly with the shock distance for both 
types of cooling function: is this due to the influence of an extended quasi adiabatic region 
where enthalpy gradients contribute to the advective-acoustic coupling (FOl) ? Or is this 
related to the cooling mechanism in the deceleration region ? How does |Q|max depend on 
the adiabatic index 7 ? Answering these questions should help us estimate the efficiency 
\Q\ in astrophysical flows with a more realistic equation of state and more elaborate cooling 
processes. 

(ii) What are the conditions for the dominance of a / = 2 mode, as observed in Fig. 6 
for a = 3/2, (3 = 5/2 for 1.5 < rsh/r^, < 1.9 ? How docs the maximum efficiency |Q|max 
depend on the degree / of the perturbation ? This question could be important with respect 
to the asymmetry of the explosion, and the subsequent kick of the neutron star (Scheck et 
al. 2006a). 

(iii) Can we better understand the conditions under which the instability is dominated 
by the fundamental mode, its first or second harmonic ? This question may be important 
with respect to the explosion mechanism proposed by Burrows et al. (2006). This mechanism 
requires the nonlinear transfer of energy from an unstable flow above the neutron star to a 
gravity mode of the neutron star, whose frequency can be a factor 10 higher. We may expect 
that the higher the frequency of the advective-acoustic cycle, the easier the excitation of the 
gravity mode of the neutron star through nonlinear processes. 

Some of these questions can be addressed by further simplifying the accretion flow in order to 
allow for analytical calculations of the advective-acoustic cycle, free from the low frequency 
limitations inherent to the WKB approximation. This is the purpose of a companion paper 
(Foghzzo et al. 2006). 
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A. Explicit relations between 5v, 5p, 5P, 5jC and /, h, 5K 

The functions /, h, 6S can be translated into the classical variables Svr, Sp, SP, Sc^ 
using Eqs. (5), (7) and (8): 



SVr 
V 

Sp 

p 

5P_ 

IP 

^2 



1 


1 - 


-M^ 




1 


1 - 


-M^ 




1 


1 - 


-M^ 


7 


- 1 


1 - 





h + 6S-l], 



f 



-M'h -SS+^j, 
-M'h-il + ij-l)M')— + ^ 



7 



- M^h - M^SS 



(Al) 
(A2) 
(A3) 
(A4) 



The functions SK, f are related to the transverse velocity perturbation (0, Svq, Sv^), accord- 
ing to the following equation, obtained from a combination of the transverse components of 
the Euler equation: 



SA 



sin 9 



d 



d 



— {sm.95ve) + 7^5f. 



de 



= -[5K-l{l + l)f]. 



(A5) 
(A6) 



Using the fact that P — pc^/j, the heating function defined by Eq. (1) and Eq. (4) is 
perturbed as follows: 



\pv 
\pvj 



- I)— + a— 

p V 



\pv J 



(A7) 
(A8) 



In these equations, the perturbations 5vz, 5p and 5c^ can be replaced by functions of /, h, 5S 
using Eqs. (A1-A4). 



B. Shock boundetry conditions 

The boundary condition at the shock follows the conservation of mass flux, momentum 
flux and energy flux in the frame of the shock: 



pi {vi - Av) = (psh + Spsh) {vsh + Svsh - Av) , 



(Bl) 
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piivi - Avf + pi — 
7 



{vi - Av f 



+ (Psh + Spsh) 



7-1 



7 

I ^ 



7 



(B2) 
(B3) 



where quantities are measured at the position + A^. Keeping the first order terms, and 
using the definition of /, h, these equations are rewritten at the position Tgh using a Taylor 
expansion: 



pivihsh - (Psh - pi)Av = A( 



7 



7 



AC 



AC 



dr 



■ 7--/1 

From the equations (2), (3), (4) of the stationary fiow, 
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We obtain: 
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(B8) 
(B9) 
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In the entropy equation, the gravity term GM/r'^^ and the term I'li'sh/'^sh due to the spherical 
geometry can be rewritten as follows: 

GM_^vjv^^ A_2^, (B13) 



Vl / Vi 



2rsh 

The transverse velocity immediately after the shock is deduced from the conservation of the 
tangential component of the velocity, in the spirit of Landau & Lifschitz (1989), leading to 
Eqs. (22-23). M^h is deduced from its definition (A5) and Eqs. (22-23): 

= -l{l + l){v,-v,^)AC. (B15) 

SKsh is deduced from SAgh using Eq. (A6) and (B12). The assumption that jCi <C >Csh, with 
Ah = Psht'shCghV-S'sh/T, leads to Eqs. (15), (16) and (17) if the shock is strong. 



C. Projection of perturbations on acoustic and advected waves 

C.l. Uniform adiabatic flow 

In a uniform, adiabatic flow moving at constant velocity in the direction z, any pertur- 
bation /, h, 5S, 5K associated with the frequency cu and perpendicular wavenumber k± can 
be decomposed as a sum of acoustic waves and advected waves as follows: 

/ = r+r+f+f", (ci) 

h = h-^ + h- + + h^, (C2) 

where an acoustic wave is noted f~^,h'^ if it propagates in the direction of the flow, and 
/~, h~ otherwise. The contribution to /, /i of a vorticity perturbation 5K such that SS — 
is: 

^ ~ 1-//W kl' 

fK 

= ^, (C4) 

V 

where /i^ = 1 — k^c'^[l — Ai'^)/uj'^. An entropy- vorticity perturbation such that 6K — 
contributes to the perturbation /, h as follows: 

fS 6S 



(? 1 - [P'M'^ 7 



(C5) 



= ^f-8S. (C6) 
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The acoustic component in a uniform adiabatic flow is deduced from f,h,5S,5K through: 



'M c 



C.2. Extension to a spherical flow with cooUng 

In a spherical flow with gradients, the advected and propagating waves are no longer 
independent, but coupled even if the flow were adiabatic. Moreover, cooling processes are 
responsible for an additional coupling between advected and propagating perturbations. We 
choose to use the same decomposition obtained in a uniform adiabatic flow, adapted to 
spherical coordinates by replacing /^^ by the spherical value (Eq. 14), and kj_ by /(/ + l)/r^ 
(the eigenvalue of the Laplacian operator) : 

^ - 1(1 + iy ^^""^ 

^ i^, (CIO) 

c2 - 1 - pL^M^ 7 ' ^ ^ 

= ^f-5S. (C12) 



. i/.|!(....)-i±fl(/^.-^). (CIS) 

For any perturbation /, h, SS, SK in a spherical nonadiabatic flow, the quantities (/=*=, f^) 
and (h^^, , h^) deflned above naturally satisfy: 

/ = ^+^+/^+/^ (C15) 

h = h+ + h- + h^ + h^. (C16) 

This decomposition describes the amount of advected and propagating waves that would be 
measured if the perturbation were allowed to continue its evolution in a uniform adiabatic 
flow. We choose to apply this decomposition at the shock radius, on the subsonic side. The 
identiflcation of with acoustic waves is strictly valid in the WKB approximation, when 
the wavelength is short compared to the scale of the gradients in the flow. The threshold of 
validity of this approximation is evaluated in Sect. 4.2.4. 



D. Numerical procedure to calculate TZy, Qv and 

One consequence of cooling is the fact that vorticity and entropy perturbations produced 
at the shock do not satisfy 6Ksh = as in adiabatic flows (Eq. 18). This raises the question 
of the acoustic feedback produced by the advection of the residual value of (5A'sh generated 
by cooling at the shock. The global efficiency associated with this feedback is thus also 
computed as a check of consistency: this feedback is indeed small in all our calculations 
(IQ^I <0.1). 



D.l. Numerical calculation of TZy, Qy, 

For a given frequency cUr and degree the differential system is integrated four times 
from Tgh to r* with the following boundary conditions: 

(i) acoustic wave propagating downward: 

SK,^ = 0, 6S,^ = 0, U = /st = 1, ^sh = (Dl) 

(ii) acoustic wave propagating upward: 

5K,^ = 0, 5S,^ = 0, U = /sh = 1, ^sh = -XT (D2) 

(iii) entropy /vorticity wave advected downward: 

5K,^ = 0, 5S,i, = 1, U = /si ^sh = hi. (D3) 

(iv) vorticity wave advected downward: 

5X,h = 1, SS,i, = 0, U = fsi, ^sh = hg. (D4) 

In each of these four cases, the differential system is integrated from the shock down to the 
accretor surface r*, where the velocity perturbation reaches a value {6v/v){r^) noted a+, 
a_, as and ax respectively. A linear combination of a couple of these integrated solutions 
allows us to construct three solutions which fulfills the boundary condition {Sv/v){r*) — 
and measure at the shock radius the following efficiencies of acoustic feedback within the 
flow: 



(i) Acoustic reflection, without any advected perturbation at the shock: 
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(ii) Acoustic feedback produced by an entropy/vorticity perturbation such that SKgh 

0: 



. -j-^^. (DO) 



(iii) Acoustic feedback produced by a vorticity perturbation such that SSsh = 0: 



D.2. Calculation of 7^sh, Qsh, 
At the shock, the couphng coefficients TZsh, Qsh and are defined by 

7^sh = ^, (D8) 

Jsh 

Qsh = ^, (D9) 

Jsh 

Ql ^ (Dio) 

/sh 

where /^j^, /J^ and /^^ are deduced from Eqs. (C3), (C5) and (C7), with the boundary values 
fsh, hsh, SSsh and 6Ksh established in Eqs. (15-18). Although cooling is neglected in the 
projection of the perturbation / on f^, (Eqs. C9-C14), some effect of cooling on the 
jump conditions is taken into account through f^, f^^ and f^. 
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Fig. 11. — Radial profile of the velocity gradient 27rdf/dr, for different stationary fiows, for 
both set of cooling parameters a = 3/2, /3 = 5/2 (thin dotted lines) and a = 6, /? = 1 
(full line). Each curve is labeled by the ratio Vsh/Ty. The velocity gradient is normalized by 
ksh|/(rsh -r*). 
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Fig. 12. — Comparison between the number riosc of eigenmodes per eigenspectrum oscillation, 
and the ratio of timescales Tv/Tac(r*), for the two cooling functions, with rgh/Tv = 30. 




Fig. 13. — Comparison between the eigenfrequencies computed from the boundary value 
problem and the first order estimate (log Q)/tv, in a flow with a — 3/2 and /3 = 5/2, for 
different shock distances Tsh/r^. The dotted lines correspond to the minimum and maximum 
growth rates described by Eq. (29), in which tq/t-jz is approximated by Ty/T^{r^). 
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Fig. 14. — Same as Fig. 13, with a different cooling function: a — 6 and P — 1. 
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Fig. 15. — Ratio of the frequencies of the first two eigenmodes a;r(l)/<^r(0) corresponding to 
I — 1 perturbations: the oscillation period li^jiOr of the fundamental mode is an excellent 
measure of the advective-acoustic timescale tq. 
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Fig. 16. — Range of unstable frequencies of the / = 1 mode. The frequency of the most 
unstable mode (thick full line) corresponds to one of the first harmonics noted "f, "hi", 
"h2" . The cut-off frequency (short dashed thick line) is a steeper function of the shock radius 
in the flow with a < (5. 
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Fig. 17. — Maximum efficiency |Q|wkb of / = 1 pcrtm^bations in the WKB approximation 
(thick hnes), compared to exp(u;jrQ) of the low frequency eigenmodes "f, "hi", and "h2". 
The cycle timescale is approximated as tq ~ 27i{k + l)/ujr{k). 
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Fig. 18. — Radial profiles of the pressure perturbation \5P\/P for the fundamental / = 1 
mode ("f) and its first harmonics ("hi", "h2"...), in the fiow with a = 3/2, p = 5/2, 
rsh/r* = 5 (upper plot) and rgh/T* = 10 (lower plot). The pressure profile is shown as a 
function of the fractional radius, with the corresponding Mach number shown on the right 
axis. In the flow with Tgh/r* = 10, the frequency of the third harmonic "h3" is high enough 
to estabhsh the instabihty of the ad vective- acoustic cycle |Q| = 4.4 and the stability of the 
purely acoustic cycle \TZ\ — 0.3 (upper plot in Fig. 8). The shapes of all these eigenfunctions 
look very similar, suggesting a common physical mechanism. A lower bound of the position 
of the effective coupling radius Ves is indicated by a circled symbol "+" , corresponding to 
Tadyii^) = 27r(A; + l)/u!r for the k*^ harmonic. 



